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ABSTRACT 

We present the results of one and three-dimensional radiative transfer calculations of polarized 
spectra emerging from snapshots of radiation magnetohydrodynamical simulations of the local vertical 
structure of black hole accretion disks. The simulations cover a wide range of physical regimes relevant 
for the high/soft state of black hole X-ray binaries. We constrain the uncertainties in theoretical 
spectral color correction factors due to the presence of magnetic support of the disk surface layers and 
strong density inhomogeneities. For the radiation dominated simulation, magnetic support increases 
the color correction factor by about ten percent, but this is largely compensated by a ten percent 
softening due to inhomogeneities. We also compute the effects of inhomogeneities and Faraday rotation 
on the resulting polarization. Magnetic fields in the simulations are just strong enough to produce 
significant Faraday depolarization near the spectral peak of the radiation field. X-ray polarimetry 
may therefore be a valuable diagnostic of accretion disk magnetic fields, being able to directly test 
simulations of magnetorotational turbulence. 

Subject headings: accretion, accretion disks — black hole physics — polarization — X-rays:binaries 



1. INTRODUCTION 

The most well understood accretion flow state onto a 
black hole is that of a geometrically thin, optically thick 
disk. Considerable theoretical effort has been devoted 
to calculating spectral models of such disks for com- 
paris o n with observations f e.g. iKolvkhalov fc Sunvaeyl 
1981 iLaor fc Netzei] fl98l iShimura &: Takaharal 1X995 



Smcell & Krolik '1998; 'Hubenv et al." 'MT; ' Davis et all 
2005; Hui et al. 2005; Davis & Hubcny 200^. These 
models have been compared to observed colors and 
spectra of ac tive galactic nuclei ( e.g. iLaod Il99d : 
iBonning et al.ll2007l:lDavis et al.ll200 



^v^..^..x.^ ■^....^wvy,,, , ultraluminous X- 

rav sour ces (iHui fc Kro hk"2008f ), and black ho le X-ray 
binaries (|Davis et al.l 12006: Don e fc Daviil2008[ ). In the 

last case, where the models appear to perform quite well 
for the high/soft state, attempts have been made to use 
them to measure the spi ns of black holes from observed 
X-ray continuum s pectra ( ghafee et al.''2006'; 'Davis et al.l 

2006; Middlcton ct al. 200l; 



ray 

20061: McClintock ' et al.l 



Liu et al.ll2008l : lGou et al 



ll009). 



It is important to bear in mind, however, that all such 
models are based on certain ad hoc assumptions. First, 
they generally adopt radial profiles of surface mass den- 
sity and emissi vity based on some form of t he alpha stress 
prescription of Shakur a fc SunyaevI |1973') and some in- 
ner boundary condition on the stress (usually zero). 
These models also make various assumptions about the 
vertical structure of the disk at each radius. In partic- 
ular, at a given radius the structure depends only on 
height above the midplane, and inhomogeneities gener- 
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ated by instabilities and turbulence are neglected. The 
disk is usually assumed to be vertically supported against 
the tidal gravity of the black hole by gas and radiation 
pressure alone. Some prescription for the vertical distri- 
bution of dissipation is adopted, e.g. that the dissipation 
rate per unit mass is constant. Finally, some mechanism 
for vertical heat transport is assumed, e.g. radiative dif- 
fusion and/or convection. 

Significant progress in understanding the vertical 
structure of accretion disks has been made with verti- 
cally stratified shear ing box simulations of magnetorota- 
tional turbulence (Balb us fc IIawlevl ll998). In particu- 
lar, radiation magnetohydrodynamical simulations of the 
vertical structure have now been done for a broad range 
of midplane radiation to gas press ure ratios of relevance 
to bl a ck hole accretion disks ([Turn er 2004; Hiros c et al.l 
2006; Krol ik et al.ll2007t iBlaeseTaLl l2007i : iHirose et all 
200 9). These simulations have demonstrated that mag- 
netic forces contribute significantly to the vertical hy- 
drostatic balance of the disk, and can be dominant 
over gas and radiation pressure gradient forces in the 
surface layers where the spectrum is formed. At the 
same time, these magnetically dominated regions are 
Parker unstable, and as a result, very large density in- 
homogeneities form near the thermalization and scat- 
tering photospheres. The mass density falls off steeply 
with height above the midplane, and most of the ac- 
cretion power is (numerically) dissipated at high opti- 
cal depth near the midplane regions. Radiative diffu- 
sion completely dominates the vertical heat transport, 
and the fraction of accretion power carried out to the 
photospheres by Poynting or mechanical fluxes is negli- 
gible. Finally, the overall turbulent stress scales approx- 
imately with the total thermal pressure near the mid- 
plane. Nonetheless, the inner radiation pressur e dom- 
inated regions of the d isk are thermally stable (|Turneil 
[200l IHirose eTall lMig^ . They may, howe ver, b e subject 
to inflow instability (Lightman & Eardle^ 19741). 

In previous work ([Davis et al.l 12005 : iBlaes et al.l 
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TABLE 1 

Parameters of Simulation Snapshots 



Simulation 


M 




Q 


Epoch 


above or below 


mo 










(rad s"-"-) 


(orbits) 


midplane 


(g cm~2) 


(K) 


0326c(=') 


6.62 


300 


5.9 


60 


below 


5.10 X 10* 


5.25 X 10^ 


0528a('') 


6.62 


150 


17 


90 


above 


4.75 X 10-* 


1.16 X 10^ 


1112a('=) 


6.62 


30 


190 


200 


below 


5.37 X 10-* 


3.45 X 10** 



IHirose et"all 1 120061) '^ iKrolik et al.l I I2007D : [slaes et al.l 1 120071) '^ IHirose et al.l I I2009D 



I2006D . we incorporated the numeri cal dissipation pro- 
files comput e d in t he simulations of iTurneil (|2004D and 
IHirose et all (|2006( ) in one dimensional non-LTE vertical 
structure models of disk annuli. For models which ne- 
glect magnetic support, we found very little difference 
in the emergent spectra compared to models based on 
the assumption of constant dissipation per unit mass, 
the common prescription which gives constant density in 
radiation dominated regions. Even though the numeri- 
cal dissipation profiles have the dissipation per unit mass 
generally increasing outward near the photosphere, most 
of the accretion power is still dissipated and thermalized 
at high optical depth, resulting in very similar spectra. 

In contrast, the addition of substantial magnetic sup- 
port of the outer layers of the disk resul ts in significan t 
differences in the emergent spectra (Blaes et al.ll2006[ ). 
Compared to mo dels that neglect magnetic forces (e.g. 
iDavis et"ani2005f ). magnetically supported atmospheres 
have larger density scale heights, resulting in smaller 
densities at the photosphere. This increases non-LTE 
effects and the ionization state of the gas, and also in- 
creases the ratio of scattering to absorption opacity. All 
of these effects result in a hardening of the dis k spectrum 
(|Blaes et al.ll2006l : iBegelman fc Pringlell2007f ). 

However, the strong density inhomogeneities that are 
seen in the surface layers may act to soften the spec- 
trum, both by increasing the thermalization of the ra- 
diation through the enhanced absorption opacity in the 
denser regions, and by enhancing the rate of cooling as 
photons di ffuse faster through the low density regions 
(jPavis et al . 2004; Bcgclnian 2006). 

The fact that accretion disks in black hole X-ray bi- 
naries in the high/soft state are expected to be elec- 
tron scattering dominated implies tha t their thermal 
spe ctra should be significantly polarized (^ Chandrasekhad 
Il96 0i). However, the complex photospheric geometry 
that results from the strong density inhomogeneities seen 
in th e simulations may reduce the degree of polariza- 
tion (jColeman fc Shield"^ Il990( l. Even if the disk sur- 
face is flat, rela tivistic effects can dilute the polarization 
(jStark &: Conn ors 1977). On the other hand, relativistic 
effects can also cause photons emitted by one part of the 
disk to scatter off a different part of the disk, enhancing 
the polarization and rotating its angle (jAgol fc KrolikI 
l200d: ISchnittman fc KrolikI [2009f ) . Strong photospheric 
magnetic fields can also Faraday depolari ze the emergent 
radiation field (iGnedin fc Silant'evlll978D . and the fields 
seen in the si mulations a.re ma rginally strong enough to 
do just that (jBlaes et al.ll2007l ). Future X -ray polarime- 
try satellites may therefore provide valuable diagnostics 
of surface magnetic fields in black hole accretion disks 
(jCnedin et all 120061 : iSilant'evI [200l . providing a good 



observational test of the simulation results. 

In this paper we present the results of spectropolari- 
metric radiative transfer calculations through representa- 
tive simulation domains that cover the range of radiation 
to gas pressure regimes that are relevant to black hole X- 
ray binaries. We quantify the degree of spectral harden- 
ing due to magnetic support, and the degree of spectral 
softening due to density inhomogeneities. We also calcu- 
late the expected polarization signatures of both of these 
effects, and confirm that X-ray polarimetry will indeed 
be able to provide interesting diagnostics on black hole 
accretion disk magnetic fields. 

This paper is organized as follows. In ^l2]we briefly de- 
scribe the numerical simulation data that we use in our 
radiative transfer calculations. In 531 we compute one- 
dimensional non-LTE atmosphere models based on hori- 
zontally averaged dissipation and magnetic force profiles 
taken from the simulations. Then in 21 we incorporate 
the fully three-dimensional structure of the simulations 
by using Monte Carlo radiative transfer calculations to 
compute the local emergent spectrum. We keep track 
of the polarization in these calculations, and we discuss 
the emergent polarization spectrum in fJS] Using certain 
scalings that we derive from comparing the different sim- 
ulations, we present an illustrative full disk spectrum of 
the polarization in ^JS] In Sj3 we discuss our results, in 
particular quantifying the uncertainties in spectral color 
correction factors and discussing how X-ray polarimetry 
could be used as a diagnostic of disk magnetic fields. We 
summarize our conclusions in iJS] Our Monte Carlo cal- 
culations fully incorporate polarized Compton scattering, 
and we summarize the relevant equations and methods in 
an Appendix. Readers interested mainly in observational 
applications of our work may wish to begin by reading 
gtland ga 

2. SIMULATION DATA 

We use data from three vertically stratified shear- 
ing box simulations in this paper: a simulation of a 
gas pressure dominated box (0326c. IHirose et al.ll200i l. 
a box with com parable midplane gas and ra,diatiq n 
pressures (;0528a. IXrolik et al.l [20071 felaes et "all 120071 ). 
and a radiation p ressure dominated simulation (1112a, 
IHirose et al.ll2009l ). The reader should consult these pa- 
pers for details of these simulations. We summarize the 
most relevant parameters for spectral modeling in Table 
1. In particular, M is the mass of the black hole, assumed 
to be non-spinning; r is the radius in the disk where the 
box would be located, in units of the gravitational radius 
rQ = GM/c^; and is the local angular velocity of the 
flow, used to approximate the local tidal gravitational 
acceleration through j^l — il'^\z\, where z is the height 
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Fig. 1. — Horizontally-averaged, vertical dissipation profiles for 
the particular epochs of the three simulations we use in this pa- 
per. In the top panel we plot the column mass density derivative 
of the vertical radiative flux, scaled with the total emergent flux 
Fo = crT^jj, as a function of column mass density m scaled with 
the total column mass density mo. The dashed, dotted, and solid 
curves correspond to simulations 0326c, 0528a, and 1112a, respec- 
tively. The bottom panel is the same as the top except that we 
plot the normalized, horizontally averaged numerical dissipation 
per unit mass e, as computed in the simulations. The dot-dashed 
lines in each panel indicate a power law with exponent 1/2, which 
corresponds to the dF/dm profile assumed in the TLUSTY calcu- 
lations. 

above the midplane. Because we are interested in the 
full three dimensional structure of the medium, we do 
not use time-averages but instead select three represen- 
tative epochs in each simulation. Table 1 also shows the 
column mass density mo to the midplane and the effec- 
tive temperature Toff of the local emergent flux at that 
epoch as computed by the simulation. 

Figure [T] shows the distribution of column mass times 
the dissipation per unit mass as a function of column 
mass m for each of the three simulation data sets. We 
measure this in two ways. For the upper panel we dif- 
ference the vertical radiative flux with respect to column 
mass, and for the lower panel we show the actual nu- 
merical dissipation rate per unit mass e. If local radia- 
tive equilibrium held exactly, then these two quantities 
should have identical distributions. This is not exactly 
true in the simulations, because a small fraction of the 
dissipated heat is transported vertically outward through 
advection of internal energy (jHirose et al.ll2009D . How- 
ever, this is most important near the midplane, and is 
completely negligible in the low column mass density re- 
gions where the emergent spectrum is formed. 

We have scaled the data shown in Figure [1] with the 
midplane column mass density mo and the emergent ra- 
diative flux _Fo for each simulation. With this normaliza- 
tion, it is remarkable how similar all the distributions are, 
despite the different physical conditions in each simula- 
tion and the different numerical resolutions in the vertical 
direction. Within the noise, all the dissipation profiles 



Fig. 2. — Horizontally-averaged, vertical magnetic accelera- 
tion amag as computed in the simulations (top) and renormalized 
magnetic acceleration aj^ag (bottom), scaled with the local tidal 
gravitational acceleration g, as a function of column mass density 
m scaled with the midplane column mass density mo. As in Fig. 
[T] the dashed, dotted, and solid curves correspond to simulations 
0326c, 0528a, and 1112a, respectively. 

are well fit by a power law of the form 



dF Ff) f m \ 
-m— — = — I I 



dm 



1/2 



mo J 



(1) 



Note that this differs slightly from the broken power 
law profile used in Blaes ct al. (2006) (Eq. 1 in that 
paper), which was based on a fit to the dissipation pro- 
file computed from the 0326c s imulation, bu t is eq uiv- 
alent to the scaling derived by iKrolik et al.l (|2007f ) for 
the 0528a simulation and used for the initial condition of 
the radiation don iinated run 1112a (Hirosc ct al. 2009). 
iBlaes et al.l ()2006D found mdF/dm oc m at high col- 
umn mass densities and oc m at low column mass den- 
sities. Since this is only an approximate relation, and 
subject to differences in averaging of the simulation pro- 
files, the small difference in the exponents is not signifi- 
cant. We adopt the 1/2 exponent because it is the best 
match to the high column density behavior in all three 
simulation, to within the uncertainties. It is important 
to recognize that all these dissipation profiles have most 
of the dissipation occurring at high column mass densi- 
ties, much deeper than the spectral forming regions of 
the annuli. 

Figure [2] shows the horizontally-averaged profiles of 
magnetic acceleration Omag, scaled with the local grav- 
itational acceleration as a function of scaled column 
mass density from the three simulation data sets. The 
magnetic acceleration is the net result of magnetic pres- 
sure and tension forces, which are comparable in mag- 
nitude in the low column mass density regions. Near 
the surface, the magnetic acceleration can be larger than 
the gravitational acceleration (even after averaging) due 
to deviations from hydrostatic equilibrium, and even be 
negative (inward). 
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Fig. 3. — Temperature (top) and density (bottom) as a function 
of height above the midplanc for an annulus. The input parameters 
have been chosen to match the t = 200 epoch of the 1112a simula- 
tion. The solid curve corresponds to a standard model, in which the 
dissipation rate per unit volume is locally proportional to density 
and magnetic support is neglected. The dashed curve also neglects 
magnetic pressure support, but incorporates the simulation-based 
dissipation profile of equation JlJ. The dot-dashed curve includes 
both the modified dissipation profile and the magnetic support. 
For each curve, the squares denote the location of the effective 
photosphere at energy 2 keV, which is near the peak in uFi, for 
the spectra shown in Fig. |4] For comparison, we also plot horizon- 
tally averaged density and temperature profiles from the simulation 
(dotted curves). 

Such excursions may give rise to interesting variabil- 
ity, but we are primarily interested in the time averaged 
behavior. In an attempt to better enforce hydrostatic 
equilibrium, we rescale the magnetic acceleration via 



Itot 



-9, 



(2) 



where atot is the total acceleration due to radiation, mag- 
netic fields, and gas pressure gradients in the simulation. 
This rescaling is also shown in Figure O It significantly 
improves the agreement between the vertical density pro- 
files of the simulation data and those produced in our one 
dimensional atmosphere models, to which we now turn. 

3. ID SPECTRAL MODELS: THE EFFECTS OF 
DISSIPATION AND MAGNETIC SUPPORT 



We use the TLUSTY code (|Hubenv fc LandflQQSh to 
compute one dimensional vertical structure models and 
emergent spectra from annuli whose input parameters 
(mo. Toff, and 17) correspond to the three snapshots sum- 
marized in Table [1] We also incorporate the vertical dis- 
sipation and magnetic support profiles measured from 
the three simulation epochs. 

The effects of the vertical support and dissipation pro- 
file are qualitatively similar in all three cases, so we only 
plot the model corresponding to 1112a in Figure [3l (A 
simi lar plot for the 03 26c snapshot can be found in Fig. 
6 of lBlaes et aLll200l ') 



The solid curve shows the "standard" annulus with 
constant dissipation rate per unit mass (i.e. constant 
dF/dm) and no magnetic support. Since the annulus 
is radiation pressure dominated throughout most of its 
extent, the density is very nearly constant, as is com- 
monly assumed. However, the effective optical depth of 
unity (marked by the squares) shows that the spectral 
forming layer occurs close to the surface where gas pres- 
sure gradients dominate and balance gravity to ensure 
hydrostatic equilib r ium. As discussed in p r evious work 
(iDavis et all l2005t iDayis fc Huben^ [20061: iDavis et all 
l2006l : lBlaes et al.ll200fil : ID one fc Daviali20(M ). this means 
that the spectral properties are sensitive to the gas pres- 
sure scale height near the surface, but rather insensitive 
to midplane properties. 




Energy (keV) 



Fig. 4. — Local emergent spectrum for an annulus viewed from 
an inclination of 55° to the surface normal. The annulus param- 
eters have been chosen to match the t = 200 snapshot of the 
1112a simulation (Pgas < Pra.d)- The curves correspond to emission 
from models with the standard dissipation and no magnetic sup- 
port (solid curve), modified dissipation but no magnetic support 
(dashed curve), and modified dissipation and magnetic support in- 
cluded (dot-dashed curve). 

A comparison of Figures O and [J] demonstrate this lack 
of sensitivity. The dashed curve corresponds to a model 
that includes the modified dissipation profile, but still 
neglects magnetic support. Figure [3] shows a significant 
change in the density profile, which is no longer constant 
at large depths and yields a higher central density. Also, 
the temperature is generally lower at a given height be- 
neath the effective photosphere relative to the standard 
model. However, the observed spectra in Figure 3] remain 
very similar with only a small amount of hardening for 
the model with enhanced surface dissipation. This can 
be understood by noting that the temperature and den- 
sity at the effective photosphere (denoted by the squares) 
both occur near the surface where the density declines 
rapidly. Because the surface pressure scale heights are 
still very similar for the two models, the spectra are rel- 
atively unchanged even though the midplane properties 
differ. This is generically the case as long as only a small 
fraction (< 10%) of the dissipatio n occurs above the ef- 
fective photosphere (see iDone fc"Davia.2008i . for further 
discussion). All three of the snapshot models obey this 
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constraint, and inspection of Figures IH [5l and [6] shows 
that modified dissipation implied by the simulations has 
only a small effect on the spectrum. 




0.1 1.0 

Energy (keV) 



Fig. 5. — Local emergent spectrum for an annulus viewed at an 
inclination of 55° to the surface normal. The annulus parameters 
have been chosen to match the t = 60 snapshot of the 0326c sim- 
ulation (Pgas > Prs.d)- The curves have the same meanings as in 
Fig. |4] Note that the dashed and solid curves are nearly identical 
in this case. 

The model including magnetic support is shown as the 
dot-dashed curve in Figure [31 The addition of mag- 
netic support greatly increases the pressure scale height 
near the surface. The resulting profiles are in qualitative 
agreement with the horizontally averaged profiles taken 
directly from the simulation (dotted curves), although 
they are not as smooth as the simulation results. Also, 
the density is not strictly monotonic decreasing as re- 
quired for stable hydrostatic equilibrium. These incon- 
sistencies result from the assumption of hydrostatic equi- 
librium in the atmosphere calculation, even though the 
simulations can have significant, non- hydrostatic dynam- 
ics near the surface. Also, the temperature differences 
arise partly due to differences in radiative equilibrium in 
the two cases. The radiation in the simulations is handled 
through flux-limited diffusion ( Levermore &: Pomranin j 
[l981) and grey opacity which only accounts for free- free 
emission and absorption. The Id calculation solves full 
radiative transfer with bound-free opacities, but neglects 
the effects of inhomogeneities in the full 3d geometry of 
the simulations. We'll expand on this point in i )4.2l 

Due to the much larger pressure scale height at 
the surface, the density (temperature) at the effective 
photosphere is m uch lower (higher) . As discussed in 
iBlaes et al.l ()2006D . this combination of higher temper- 
ature and lower density leads to harder spectra, both 
due to the effects on absorption features at lower tem- 
peratures and due to the effects of electron scattering at 
higher temperatures. Indeed, Figures [H El and |6l show 
that magnetically supported annuli all give harder spec- 
tra (dot-dashed curves), demonstrating that the quali- 
tative result (harder spectra) is quite general. We can 
make this more quantitative by comparing the color cor- 
rected blackbody models that best match the magnetized 
and standard atmospheres. This estimate is most ro- 



bust for the model corresponding to the 1112a snapshot 
since the strong absorption edges in the other models 
prevent the color-corrected blackbody from providing a 
good fit, making the color correction (or hardening fac- 
tor) somewhat ambiguous. For the 1112a model, which 
lacks significant edge features, we find the color correc- 
tion increases from about 1.67 to 1.84 when incorporating 
magnetic support, a roughly 10% increase. 
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Fig. 6. — Local emergent spectrum for an annulus viewed from 
an inclination of 55° to the surface normal. The annulus parame- 
ters have been chosen to match the i = 90 snapshot of the 0528a 
simulation (Pgas ~ frad)- The curves have the same meanings as 
in Fig. g] 

4. 3D SPECTRAL MODELS: THE EFFECTS OF 
INHOMOGENEITIES AND MAGNETIC FIELDS 

4.1. Monte Carlo Radiative Transfer 

Although the hydrostatic. Id calculations described in 
121 capture much of the most relevant physics, they ne- 
glect all the 3d, dynamical information available from the 
simulation results. In order to explore these effects we 
have computed fully 3d radiative transfer calculations us- 
ing Monte Carlo (MC) methods for each snapshot listed 
in Table dl 

The MC code reads in a grid of densities, temperatures, 
and magnetic field vectors from a single 3d snapshot of a 
simulation. These variables are taken to be constant in 
time throughout the MC calculation. At the beginning of 
the calculation, the electron density, free-free emissivity 
and absorptivity are computed for each zone. We assume 
the plasma is completely ionized with a ten percent num- 
ber abundance of helium. (The MC calculations do not 
account for partial ionization, non-LTE effects or bound- 
free processes. These are included in the Id calculations 

of El) 

The code then propagates photon packets through the 
domain until they escape or are absorbed. The domain is 
assumed to be periodic in both the radial and azimuthal 
directions^, so escape only occurs through the vertical 

^ The simulations themselves assume shearing periodic boundary 
conditions in the radial direction (Hawlcv ct al. 1995). Because 
we neglect time evolution, we simply assume periodic boundary 
conditions in the radial direction for the Monte Carlo calculation. 
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Fig. 7. — Local emergent intensity (top) and polarization (bot- 
tom) from the t = 60 snapshot of the 0326c simulation (Pgas > 
^rad)i viewed at an inclination of 79° to the surface normal. In 
the top panel, the spectra are generated from fully 3d (black curve) 
and horizontally averaged Id (red curve) Monte Carlo calculations. 
In the bottom panel, we compare the polarization from Id calcula- 
tions (red curve), and 3d calculations with (blue curve) and without 
(black curve) the effects of Faraday rotation. The dotted curve in 
the bottom panel corresponds to fitting function given by eq. \5\ 

boundary. For efficiency, we only include the outer por- 
tions of the full disk volume near the surface, and assume 
reflecting boundaries at some inner base. We choose this 
inner boundary to be at a large optical depth to true ab- 
sorption (> 10 for typical energies) so that its location 
has negligible effect on the output spectra and polariza- 
tion. 

The initial location of each packet is chosen randomly 
throughout the domain, and each packet is assigned a 
weight based on the emissivity of the initial grid zone. 
The initial direction vector of the packet is randomly 
chosen from an isotropic distribution. The photons are 
initially unpolarized and the energy is randomly chosen 
and assigned a weight to approximate the frequency de- 
pendence of the free-free emissivity. 

After initialization, the absorption and electron scat- 
tering mean free paths are calculated for the appropriate 
energy and the packet is advanced along a ray pointing 
along the direction vector by an optical depth chosen ran- 
domly from an exponential distribution. If this optical 
depth is sufficiently large, the packet may be propagated 
into a neighboring grid zone. The code assumes each grid 
zone is homogeneous in density, temperature, and mag- 
netic field strength. The mean free paths are updated 
each time the packet enters a new grid zone, and prop- 
agation continues until the packet moves the prescribed 
optical depth or exits the domain. If the packet escapes, 
its weighted contribution is added to the output arrays 
corresponding to its final direction and energy, and a new 
packet is initialized. 

If the packet remains in the domain, its weight is re- 
duced by the ratio of absorption to total (absorption 



and electron scattering) opacity. If the packet weight 
falls below a prescribed minimum weight, it is consid- 
ered absorbed and a new photon packet is initialized and 
propagated. (The minimum weight is chosen to be suf- 
ficiently low that the absorbed photons have negligible 
effect on the output spectrum.) Otherwise, the packet is 
assumed to have scattered and the energy, direction vec- 
tor, and polarization vector are updated. We model the 
scattering process using t he MC methods for Com pton 
scattering as described in 'Pozdniak ov et al.l (1X983') , but 
with modifications to include polarization as described 
in the appendix. A new optical depth is then drawn and 
the process is repeated until the packet escapes or is ab- 
sorbed. 
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Fig. 8. — Local emergent intensity (top) and polarization 
(bottom) from the t = 90 snapshot of the 0528a simulation 
(^gas ~ ^'rad)i viewed at an inclination of 79° to the surface nor- 
mal. The dotted and binned curves have the same meaning as in 

Fig.m 

In order to test our MC routines we ran a number 
of comparison calculations with other methods. Since 
we did not have a single comparison method that could 
suitably handle all the complexity of the MC calcula- 
tion, we separated our test into two parts. We first 
tested the initialization and photon propagation routines 
by comparing with Id Feautrier calculations that model 
polarized Thomson scattering and allow for temperature 
and density variations in a single dimension. We then 
ran our MC code with polarized Thomson scattering 
on fully 3d grids with equivalent density and tempera- 
ture variations in the vertical dimension, finding good 
agreement for both the specific intensity and polariza- 
tion spectra. We next tested the polarized Compton 
scatte ring routi nes by comparing with results from a 
code ()Hsu fc Blacs 199^ based on th e itera tive scatter- 
ing method of iPoutanen fc Svensso i (fT99l) . We used 
both methods to calculate the emergent specific inten- 
sity and polarization from a homogeneous, finite optical 
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depth "corona" with a blackbody seed photon distribu- 
tion at its base, finding good agreement. 
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Fig. 9. — Local emergent intensity (top) and polarization 
(bottom) from the t = 200 snapshot of the 1112a simulation 
(Pgas < Piad)t viewed at an inclination of 79° to the surface nor- 
mal. The curves have the same meaning as in Fig. [7] 

4.2. The specific intensity: effects of inhomogeneities 

Figures [71 O and [9] show the resuhs of the MC calcu- 
lations described above for the 0326c, 0528a, and 1112a 
snapshots, respectively. For the moment, we concentrate 
on the top panels of each figure, which show the specific 
intensity at an inclination of 79°. As one would expect 
on physical grounds, the specific intensity shows little 
azimuthal variation, so we plot the azimuthally-averaged 
output of the MC calculation to improve the statistics. 
The solid black curves correspond to the results of MC 
calculations through the full 3d grid and the red curves 
are the MC spectra from horizontally averaged Id grids. 
In each case both sets of calculations yield rather similar 
spectral shapes, and the spectral peaks occur at nearly 
identical energies. The main difference is that the 3d 
model intensity is greater at all but the lowest energies. 
Although we have only shown the spectra at 79°, this re- 
sult qualitatively holds at all inclinations, implying that 
a larger flux is being radiated in the 3d calcula tions, in 
good agreement with previous results l| Davis et al.H2004h . 

This flux enhancement is primarily the result of the 
density inhomogeneities. Although there are some hori- 
zontal variations in the temperature, they are generally 
smaller than the variations in the density. The weaker 
variation in temperature results from the relatively short 
photon diffusion time scale, which ef ficiently smooths ou t 
temperature fluctuations (see e.g. [Turner et al.l l2005f ). 
The large density inhomogeneities are a result of highly 
variable magnetic forces that dominate the surface lay- 
ers. Densities range over factors of roughly ten above 
and below the mean at the horizontally averaged ef- 
fective photospheres in all the simulations (|Blaes et al.l 



l2007HHirose et al.|[2009f) . It is possible that photon bub- 
bles, if they exist, may make th e density even more inho- 
mogeneous ([Turner et al.ll2005l ) in the radiation pressure 
dominated regime. 

The density inhomogeneities enhance the total emis- 
sion because free-free emission (and thermal emission 
processes in general) scale as the square of the density. 
Therefore, even though average density is identical in the 
two calculations, the excess emission from high density 
regions more than compensates for the deficit in low den- 
sity regions. Of course, free-free absorption also scales as 
the square of the density and if it were the only opacity 
source, increased absorption would offset the increased 
emissivity. However, since electron scattering is propor- 
tional to one power of density and dominates the opac- 
ity at most frequencies, photons emitted in high density 
regions scatter out to lower densities before they are ab- 
sorbed, increasing their likelihood of escape. 

Ultimately, we would like to address whether or not 
these inhomogeneities are essential to obtaining the cor- 
rect specific intensity. The fact that both the Id averaged 
domain and 3d domain yield the same spectral shape 
suggests that the spectra may be insensitive to the in- 
homogeneities. However, the flux enhancements suggest 
we might be able to get an equivalent flux with a lower 
surface temperature, possibly lowering the spectral hard- 
ening. The problem with comparing the Id horizontally 
averaged and full 3d MC calculations is that they don't 
represent the same radiative equilibrium. In a real disk 
the flux of radiation (and therefore the radiative equi- 
librium) is flxed by the overall amount and location of 
turbulent energy dissipation. One needs to compare the 
3d calculation against a Id calculation with the same 
radiative equilibrium. Fortunately, we have already de- 
scribed such calculations in iJ3| 
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Fig. 10. — Local emergent intensity at 55° from the t = 200 
snapshot of the 1112a simulation (Pgas < Prad)^ viewed at an in- 
clination of 79° to the surface normal. The binned curve is identical 
to the black binned curve in Fig. [9] The solid, unbinned curve is a 
TLUSTY calculation similar to the dot-dashed curve in Fig. [J] but 
without bound- free opacity and with a lower effective temperature 
as discussed in i|4.2l 

Figure [TUl shows a comparison of the 3d MC spectrum 
(bins) from Figure [9l with an equivalent TLUSTY at- 
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mosphere spectrum (solid curve). The TLUSTY model 
includes the simulation derived dissipation profile and 
magnetic support, but differs from the one shown in Fig- 
ure [4] (the dot-dashed curve) in that we have dropped 
the bound-free opacity and used a slightly lower effec- 
tive temperature to better match the 3d MC calculation. 
The lower effective temperature is needed because the to- 
tal flux in the MC calculation is slightly lower than that 
found directly by the simulation. We find that the 3d 
MC spectrum is significantly softer than the Id TLUSTY 
spectrum. In fact, the peak is lower by about 12%, al- 
most completely offsetting the hardening due to the mag- 
netic support discussed in ij3l 

In fact, a similar cancellation also occurs in the 0326c 
and 0528a snapshots as well. Again, we repeated the Id 
TLUSTY calculations from |21 but neglecting bound-free 
opacities and fixing the effective temperature to match 
the MC results. The spectral hardening due to mag- 
netic support was about the same level as the softening 
due to inhomogeneities in both cases: ~ 10% for 0326c 
and < 5% for 0528a. Since the magnetic fields providing 
the support are also the primary source of the inhomo- 
geneities, a correlation between the effects is not surpris- 
ing. However, it's remarkable that the magnitude of the 
effects are such that they cancel out in all three sets of 
calculations. 
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Fig. 11. — Temperature as a function of height above the mid- 
plane for an annulus. The sohd curve is a Id atmosphere calcu- 
lation with input parameters chosen to match the t = 200 epoch 
of the 1112a simulation, but with a shghtly lower effective tem- 
perature and neglecting bound-free opacity in order to match the 
3d MC calculation. For comparison, we also plot the horizontally 
averaged temperature profile from the simulation (dotted curve). 

We plot a comparison of the temperature profiles from 
the two calculations in Figure [TlJ The dotted curve 
shows the horizontally averaged temperature from the 
1112a simulation snapshot and the solid curve is the Id 
atmosphere model. Overall, the Id profile is a better 
match to the simulation than the equivalent model in 
Figure [3] (dot-dashed curve) since we have modified the 
input parameters to better match the MC calculation. 
Nevertheless, the temperature gradients differ near the 
surface and the effective photosphere (square symbol) is 
~ 15% higher than the simulation at the same height 



above the midplane. We generated similar plots for the 
other two snapshots and they are qualitatively consistent 
with Figure [H] 

This result, along with the flux enhancements dis- 
cussed above, strongly suggests that the inhomogeneities 
allow the same flux to be radiated with average lower 
surface temperatures, leading to a softer spectrum. How- 
ever, one must note that a number of other differences 
between the two calculations may also contribute to the 
resulting discrepancies in the temperature profiles and 
spectra. The treatment of radiation transport in the 
simulations (grey opacity, flux limited diffusion) differs 
from both the MC and the TLUSTY calculations. Also, 
non-hydrostatic motions are significant in the surface lay- 
ers of the simulations while the Id TLUSTY calculation 
strictly enforces hydrostatic equilibrium. Distinguishing 
between these effects is difficult, but it is clear that the 
3d and dynamic nature of the simulations plays a non- 
negligible role in determining the disk spectrum. Ignor- 
ing these effects will lead to an underestimate of the spec- 
tral hardening. 

5. POLARIZATION: THE EFFECTS OF COMPTON 
SCATTERING, INHOMOGENEITIES, AND MAGNETIC 
FIELDS 

For viewing angles inclined to the surface normal, ho- 
mogeneous atmospheres will appear moderately linearly 
polarized (up to 12%) at photon energies for which elec- 
tron scatte ring is nearly elastic a nd dominates the opac- 
ity (see e.g. lChandrasekhaill960( l. Since electron scatter- 
ing opacity is typically dominant for the energy ranges of 
interest, the MC calculations do produce significant lin- 
ear polarization at moderate to high inclinations. This 
can be seen in Figures [7] - IHl which plot the polarization 
viewed from an inclination of 79°. 

These polarization results are most easily discussed in 
terms of the Stokes parameters corresponding to the total 
specific intensity / and the linearly polarized intensities 
Q and U. (Electron scattering does not impart circular 
polarization so the Stokes parameter V is identically zero 
in our calculations.) We will also find it convenient to 
define the total polarization P = (Q^ + J7^)^/^. In this 
paper we always plot the normalized Stokes parameters 
{q = Q/I, u = U/T) and degree of polarization p = 
(g2+M2)i/2. 

We define Q such that positive and negative Q cor- 
respond (respectively) to polarization perpendicular and 
parallel to the surface normal. U corresponds to polar- 
ization defined relative to axes rotated by 45°. For a 
homogeneous domain, U would also be identically zero 
by symmetry. Although the simulation domains are not 
homogeneous, they are still highly symmetric, and U is 
nearly zero in all three snapshots for all viewing angles. 
Therefore, we only plot Q in Figures [7] - [9] Like the in- 
tensity, we do not observe significant azimuthal variation 
in the polarization, so we plot only the azimuthal aver- 
age. The only exception is when the effects of Faraday 
rotation are included, as discussed in 

The direction of the polarization is generally deter- 
mined by the angular dependence of the emergent radia- 
tion field. Limb darkened radiation yields polarization in 
the plane of the disk (positive Q), isotropic radiation is 
unpolarized, and limb-brightened radiation provides po- 
larization parallel to the surface normal (negative Q). 
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Since Thomson scattering dominated atmospheres are 
hmb darken ed, one nominally exp ects polarization with 
positive Q (jChandrasekhail [l960[ ) . However, a number 
of other processes are relevant: absorption at low ener- 
gies, changes in photon energy due to inelastic electron 
scattering, Faraday rotation by the magnetic fields, and 
geometric effects due to the inhomogeneities. 

Absorption opacity changes the polarization both by 
competing with scattering opacity and by altering the 
anisotropy of the radiation field. Depending on the ver- 
tical gradient of the thermal source function, this can 
increase the polarization compared to a pure scattering 
atmosphere (lHarringtoiJl969tlLoskutov fc Sobolevlll979l : 
iBochkarev et al.lll985[ ). or~^crease it and even turn it 
negative (the "Nagirner effect" , c.f. iGnedin &: Silant'evI 
[197^ . For the free- free opacity assumed here, absorption 
dominates at lower photon energies, and the Nagirner ef- 
fect is clearly evident in the bottom panel of Figure [7] at 
energies below 0.03 keV. It is also present in the other 
two snapshots, but at lower energies than those plotted 
in Figures [5] and [5] 
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Fig. 12. — Local emergent intensity (top) and polarization (bot- 
tom) from the t = 90 snapshot of the 0528a simulation (Pgas ^ 
^rad)i viewed at an inclination of 79° to the surface normal. The 
black and red curves correspond to full 3d MC calculations where 
electron scattering is treated in the Compton and Thomson limits, 
respectively. The dashed curve shows the prediction for a semi- 
infinite, scattering-dominated atmosphere viewed from the same 
inclination. 



5.1. Compton scattering 

The effects of inelastic (Compton) electron scattering 
manifest themselves as a turnover in Q at photon ener- 
gies well above the spectral peak. If electron scattering 
is treated in the elastic (Thomson) limit the polarization 
would remain positive. We have deliberately neglected 
to plot this turnover in Figures [7| -[3] due to the poor 
statistics at these energies. But, the effect can be seen in 
Figure [T^ where we compare a full 3d model where elec- 
tron scattering is treated in the Compton limit (black 



curve) with one in which the it's computed in the Thom- 
son limit (red curve). 

At high energies above or near the peak, the predic- 
tions for polarization differ significantly. The result is not 
due to the differences in the scattering matrices, since the 
full Compton-scattering matrix is well approximated by 
the Rayleigh matrix at the energies of interest. (See the 
appendix for further details.) Instead, it is related to the 
different frequency dependence of the angular distribu- 
tion of the radiation in the two cases. Since the photon 
energies are fixed in the Thomson model, high energy 
photons tend to be preferentially emitted deeper in the 
atmosphere where temperatures are greater. Since pho- 
tons are initially emitted at large scattering depths, the 
radiation field becomes limb darkened, the spectrum be- 
comes a modified blackbody (Rybicki fc Lig htman 1979), 
and the polarization is well approximated by the pre- 
dictions for a homogeneous, semi-infinite, scatteri ng- 
dominated atmosphere (see e.g. lChandrasekhailll960D . 

However, when Compton effects are included these 
high energy photons emitted at large depth can now 
efficiently exchange energy with the electrons near the 
surface. Since the surface is cooler, these electrons are 
typically lower in energy and photons will tend (on aver- 
age) to lose energy until they have a mean energy similar 
to the photospheric electrons. This leads to the Wien tail 
at high photon energies seen in the specific intensity for 
the Compton spectra. The photon angular distribution 
and the energy distribution are now coupled by the scat- 
tering history, and the polarization can therefore differ 
from that of an optically thick Thomson scattering at- 
mosphere. Indeed, the polarization of 2 keV photons is 
slightly enhanced by Compton scattering in Figure 1121 

Near the surface photon diffusion smooths out most, 
but not all, temperature inhomogeneities. In a few lo- 
calized regions the gas temperature significantly exceeds 
the effective temperature of the radiation, though this 
is likely artificial as simulation 0528a neglected Comp- 
ton energy exchan ge between the radiation and the gas 
(jBlaes et al.l 120071 ). In these pockets, photons are up- 
scattered to higher energies, producing the break seen 
at high energy in the top panel of Figure [121 Due to 
the small optical depth of the hot material, these pho- 
tons only overcome the emission in the Wien tail at high 
energies. At these energies the up-scattered emission is 
modestly limb-brightened, due to the slightly increased 
probability of laterally propagating photons to be scat- 
tered before escape. This is completely analogous to the 
polarization flip produced by inverse Compton scattering 
in an optically thin, hot corona in plane para llel geom- 
etry (jHaardt fc Mattlll993t iHsu fc Blaesl |l998D, and is a 
reminder that such coronal geometries, if they exist in 
nature, will produce high energy radiation that is polar- 
ized parallel to the projected disk rotation axis. 

5.2. Inhomogeneities 

The effects of inhomogeneities can be seen in the bot- 
tom panels of Figures [7] - [51 The black curves show the 
polarization of the full 3d calculations (excluding the ef- 
fects of Faraday rotation) while the red curves correspond 
to the Id, horizontally averaged domain. We find that 
the 3d polarization curve is generally less polarized than 
the Id horizontally-averaged curve, although the effect is 
rather small. 
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Fig. 13. — Contours of constant escape fraction overplotted on 
a two dimensional density slice in the x — z plane of the 0528a 
snapshot. The white contours correspond to the fraction of upward 
moving photons that escape from a grid zone after any scattering 
or emission of a packet. From top to bottom, they correspond 
to escape fractions / = exp(— t) with r = 0.5, 1, 2, and 4. The 
plotted region corresponds to the uppermost 31% of the simulation 
domain. 

This slight reduction in polarizati on is in quaUtative 
agreem ent with the expectations of iColeman fc Shieldsl 
(|1990f ). who argued that deviations from a planar photo- 
sphere (e.g. a corrugated surface) might produce a more 
isotropic, or even limb brightened, radiation field. They 
postulated that such a photospheric geometry could, in 
principle, explain the apparent results that AGN are at 
most weakly polarized (< 2%) and that the polarization 
is perpendicular to the inferred disk plane, contradicting 
the standard expectations of scattering dominated atmo- 
spheres. However, the effect observe d in our calculations 
is mu ch smaller than suggested by IColeman fc Shieldi 
(2993) and may partly be the result of the increased ra- 
tio of absorptio n to scattering opac i ty at some energies. 

The model of 'Colem an fc Shieldsl ()199GD implicitly as- 
sumes the photosphere geometry is complex, but that 
characteristic length scales of the horizontal variations is 
significantly greater than the photon mean free path. (If 
the horizontal length scale is smaller or comparable to 
a mean free path, their assumption of a IChan drasekhail 
([1960) limb-darkening profile relative to the local surface 
normal would be invalid.) 

To compare with their model we examined the shape 
of the photosphere in our 3d calculations. In Id at- 
mospheres the photosphere corresponds to a surface of 
constant escape fraction, a concept which is easily gen- 
eralized to 3d. Therefore, we computed in each grid 
zone the fraction of photon packets escaping after each 
scattering or emission event. In Figure [13] we plot the 
contours of constant escape fraction as white curves on 
top of the density in a 2d slice of the 0528a snapshot. 
Although deviations from a planar surface are signifi- 
cant, the length of horizontal variations are comparable 
to or smaller than a typical mean free path to electron 
scattering. Therefo r e, the geometric effect discussed by 
IColeman fc Shieldi ()1990f ) is nearly negligible for these 
simulation domains. 

Nevertheless, th ese results do not invalidate the 

IColeman fc Shieldi (|1990f ) hypothesis for real systems. 



Since periodic boundary conditions are assumed, larger 
length scale horizontal variations are not possible in the 
current simulations. The surface inhomog eneities are pr i- 
marily the result of Parker instability (Bla es et al.ll2007[) . 
so larger wavelength variations may be possible in actual 
accretion flows, and this hypothesis should be reexam- 
ined when larger domains become computationally fea- 
sible. 

5.3. Magnetic fields 

The bottom panels of Figures [7] -[9] also show the effects 
of magnetic fields on the polarization. Both the blue and 
black curves are 3d calculations, but with and without 
(respectively) the effects of Faraday rotation included. 
Each time a photon packet is propagated through a 
Thomson optical depth tt , the polarization vector is ro- 
tated by an amount 

XP = |^B.k (3) 

where B is the magnetic field in the current grid zone, 
k is a unit vector in the direction of photon propaga- 
tion, e is the electron charge, and A is the wavelength. 
For the observed polarization, the rotation that occurs 
after the last scattering is typically most important, so 
r ^ 1. Since k is a unit vector, the typical rotation 
angle only depends on A, B, and the direction. Due to 
the A^ dependence in equation ([3|), one is normally not 
concerned with Faraday rotation in X-ray sources. How- 
ever, the near equipartition strengths of magnetic fields 
in the simulation (~ 10^ — 10^ G) are large enough that 
Faraday rotation can be significant even at these short 
wavelengths. 

For any single photon, Faraday rotation can only 
change the direction of the polarization, not its mag- 
nitude. However, due to departures from uniformity in 
the magnetic field and differences in trajectories, differ- 
ent photons will experience rotation through different 
angles. This is particularly significant for low energies 
where xp 3> 1 and small differences along nearby tra- 
jectories can lead to rather different polarization angles. 
As a result, photons are imparted with a nearly uniform 
distribution in polarization, yielding zero net polariza- 
tion after integration. Therefore, Faraday rotation tends 
to completely depolarize at low energies, gives modest 
reductions at intermediate energies where Xf ^ 1, and 
has negligible effect for higher energies where xf ^ 1- 

For photon energies with xf ~ 1, there is an az- 
imuthally dependent rotation, due to the nearly toroidal 
net magnetic field in the snapshot. Figure 1141 shows the 
polarization angle 

= 0.5tan-i(C//Q) (4) 

plotted as a function of azimuthal angle 4> for the 0528a 
snapshot. A "i/; of zero corresponds to polarization par- 
allel to the simulation domain surface. This plot shows 
the polarization angle for an inclination of 79°, averaged 
over photon energies from 0.3-1 kcV, where xf ~ 1- For 
(f) ^ 110° and 290° the escaping photons are traveling 
parallel and anti-parallel to the net toroidal field, obtain- 
ing non-zero polarization angles. When averaged over all 
azimuth, ip is consistent with zero. 

The observed spectrum from an unresolved disk is a 
weighted average over all azimuth. It is "weighted" be- 
cause relativistic beaming can enhance emission from the 
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Fig. 14. — Polarization angle as a function of azimuthal angle 
for the t = 90 snapshot of the 0528a simulation. The polarization 
angle is computed for an inclination of 79°, and averaged over 
photon energies from 0.3-1 keV, where the Faraday rotation angle 
is near unity. 

approaching side and diminish emission from the reced- 
ing side of a disk annulus. Therefore, cancellation will 
not be complete deep in the relativistic potential well 
of a black hole, so some net rotation may be observed. 
Nevertheless, we have chosen to plot the uniformly az- 
imuthally averaged q in Figures [7] -[H since the details of 
the relativistic beaming are complicated to model, and 
depend precisely upon where the simulation would be lo- 
cated in a real accretion flow. For higher or lower photon 
energies, there is very little net rotation, and U remains 
consistent with zero, even before azimuthally averaging 
the output. 

In the bottom panels of Figures [7] - [9] the dashed curve 
corresponds to the simple analytic relation 

g = 9o(M)[l + (A/AB)2]-\ (5) 

where 90 (m) is the non-magnetized polarization value 
appropriate for the observ ed inclination in a s catter- 
ing dominated atmosphere (jChandrasekhaiiri96Cll ). = 
87r'^e/(3i3ph), Bph is the average magnitude of B at the 
photosphere of the snapshot, and is the cosine of the 
inclination angle. The quantity (A/As)^ is roughly the 
Faraday rotation angle for unit optical depth. For ener- 
gies near or below the peak, this functional form provides 
a rather good approximation to the 3d MC results when 
Faraday rotation is included. It breaks down at high 
energies because it docs not account for the Conipton 
scattering turnover, but is a nearly perfect match to the 
Thomson scattering calculations. 

Figures [7] - [9] are all plotted for the same inclination 
of 79°, a relatively high value at the edge of the ob- 
ser ved distribution for black ho le X-ray binaries (see 
e.g. lNaravan fc McChnt"ocg i2GG5V Since polarization in- 
creases with inclination, it will generally be lower in most 
sources and the effects discussed here will be more dif- 
ficult to measure. The inclination dependence of Q is 
shown in Figure [TSl The black, blue, and dashed curves 
have the same meaning as in the bottom panels of Fig- 
ures [7] -[3 

These results can be compared with previous work 



on Id atmosphere s with uniform (pilant'evi 120021 : 
IShternin et "al]l2003t iGnedin et all 120061 ) and~turbulent 
(|Silant'evll2007D rnagnetic fields. Our ad hoc relation (O 
is qualitatively similar to the polarization dependence de- 
rived in these works. They find depolarization at photon 
energies with large Faraday rotation angles and make de- 
tailed predictions for the angular and energy dependence 
of the polarization. Their analytic formulas are most 
precise in the asymptotic limit of large Faraday rotation 
angles, but it is difficult for us to probe the asymptotic 
dependence due to limited photon statistics at low ener- 
gies. For uniform magnetic field, they find p oc (5~^ for 
5 1, where 5 is half the Faraday rotation angle for op- 
tical depth of unity. This is the same depende nce as in 
( 5l) fo r A ^ As. For isotropic turbulent fields. [Silant'evI 
( 2007i ) predicts p oc A"'* for large A, a somewhat stronger 
dependence than found here. Although the simulation 
magnetic fields do have turbulent fiuctuations, there is 
a mean toroidal field which may explain our consistency 
with their uniform field relation. 
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Fig. 15. — Polarization viewed from four different inclinations 
to the surface normal. The spectra are generated from fully 3d 
calculations using the 4 = 90 snapshot of the 0528a simulation 
(Pgas ^ ^Vad). The black and blue curves represent calculations 
that exclude and include the effects of Faraday rotation, respec- 
tively. The dotted curve corresponds to the fitting function given 
by eq. [5] 



6. FULL DISK POLARIZATION MODEL 

The effects discussed in Sj5]indicate that the assumption 
of I Chandrasekhar (1960) polarization for a semi-infinite 
atmosphere is not correct at either high or low photon 
energies. However, it is important to keep in mind that 
our simulations only represent a small patch of the disk, 
and (at best) only approximate the emission from a sin- 
gle, narrow annulus in the accretion flow. In order to 
assess the possibilities for obtaining constraints from ac- 
tual observations, we need to compute full disk models. 
Given the uncertainties and complexities involved, our 
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goal is only an approximate model, which is presented to 
motivate future calculation s and observations. 

Following previous work (jHubenv et al.|[2000l and ref- 
erences therein) , we assumed a fully -r elativistic thin disk 
model (|Shakura fc SunvaevI 119731 : iNovikov fc Thornd 
[1973) with zero torque inner boundary condition at the 
innermost stabl e circular o rbit. We also use the KER- 
RTRANS code (jAgoll Il997l ) to compute the relativistic 
effects on photon geodesies and parallel transport the 
polarization vectors as they travel from the disk surface 
to the observer at infinity. 

For the sake of simplicity, we will assume that the 
local spectra are color-corrected blackbodies with a 
color-correction of 1.7, which is a qualitatively rea- 
sonable approximation at the accretion rates of inter- 
est (ShimuraJ^TaJtaims liW5; J3^^ Of 
course, the results presented in §3] (and previous work) 
indicate that this will not be quantitatively correct in de- 
tail. The main obstacle to a more precise calculation is 
that we do not have simulations for the innermost rings 
which account for most of the radiation, so we can not 
robustly model the effects of dissipation and magnetic 
support for the whole disk. Based on the simulations 
presented in this paper (which, however, may not gener- 
alize to hotter annuli), the overall shift in the spectrum 
may be very small, and the numerous approximations 
(outlined below) ultimately make this a qualitative en- 
terprise, anyway. 

We additionally assume that the local spectra have 
an angular distribution that approximately matches the 
[Chandrasckhar (1960) limb-darkening profile, although 
this only has a small effect on the output spectrum. Note 
that we are also assuming that the local rest frame po- 
larization is azimuthally symmetric and either parallel 
or perpendicular to the disk plane. We neglect the az- 
imuthal dependence of polarization position angle dis- 
cussed in section 5.3 above, due to the much larger un- 
certainties inherent in this illustrative full disk model. 

For the energy dependence of the polarization, we use 
equation ^ as an analytic approximation for each indi- 
vidual annulus. This accounts for the effects of Faraday 
rotation, but omits two other polarization-changing ef- 
fects that could be important. If there is a hot, optically- 
thin corona above the disk surface, Compton scattering 
within the corona could rotate the polarization angle to 
be nearly perpendic ular to the surface at energies above 
the t hermal peak (jHaardt fc MattI Il993t iHsu fc Blaei 
Il998| ) . In the simulation data we study, the gas in the re- 
gion above the thermalization photosphere is hotter than 
the effective temperature, but the inverse Compton scat- 
tering occurring there is much too weak to explain the 
coronal emission typically observed; consequently, our 
calculations likely underestimate the importance of this 
effect, i.e., overestimate the energy at which the polar- 
ization rotation takes place. Secondly, we do not con- 
sider the scattering of returning ra diation (jAgol fc KrolikI 
I2000t ISchnittman fc KrolikI [2009( 1. a global general rel- 
ativistic effect that also rotates the polarization (and 
strengthens it) at energies above the maximum effective 
temperature found in the disk. 

The last requirement is a specification for _Bph as a 
function of radius in the disk model. In the three 
shearing-box simulations, i?ph correlates well with the 
equipartition magnetic field strength Bcq, the magnetic 



field strength for which magnetic pressure equals the to- 
tal midplane pressure. Although there is some scatter in 
the correlation, we find that i3ph — i?Gq/40. Therefore, 
we estimate Bph by first computing B^q as a function of 
radius in our thin disk models and then use this relation 
to obtain Sph. 

This estimate is particularly uncertain, principally be- 
cause the total midplane pressure depends on a in the 
thin disk model, so that our estimate of i?ph inherits 
this dependence: Bph oc a~^/^. We assume a = 0.01 
in our models, a value which is ~ 1/2 the typical long- 
term time-average found in our stratified shearing box 
simulations, and roughly an order of magnitude smaller 
than both the typical numbers found in the disk body in 
global disk sim ulations and som e observationally-based 
estimates (e.g. iKing et all 120071 ). Using a constant a 
for this purpose may create a further problem in that 
global disk simulations often show a sizable increase in 
the time-average of this quantity in the region just out- 
side the ISCO, where a significant part of the luminosity 
is created. An additional uncertainty is introduced by 
the fact that what matters for Faraday rotation is B • k, 
and the relative magnitudes of the different components 
of B may change systematically with radius. 
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Fig. 16. — Full disk specific intensity (top panel), degree of po- 
larization (middle panel) and polarization position angle (bottom 
panel), as viewed by an observer at infinity at an inclination of 73° 
to the disk normal. All spectra are generated from relativistic disk 
models assuming L/L-^dd = 0.1, a = 0.01, and M = IOMq, which 
are described further in the text. The solid and dashed curves are 
computed for a, = and a* = 0.99, respectively. In the bottom 
panels, the thick curves correspond to our fiducial model for which 
eq. [5] is used to specify the polarization emitted from each annu- 
lus and Bpi^ = iJcq/40 is assumed. The thin curves correspond to 
models where the polarization is assumed be the Chandrasckhar 
value. 

After combining these prescriptions, we compute full 
disk models, the results of which are shown in Figures fTHl 
andfTTl The top panels show the total specific intensity, 
the middle panels show the degree of polarization, and 
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the bottom panels show the polarization angle tp. Po- 
larization parallel to the plane of the disk corresponds 
to ip = 0°, while ip = 90° corresponds to polarization 
parallel to the projected rotation axis of the disk. 



It has been well established by previous work that the 
polarization is sensitive to the pr operties of the spacetime 
and observer inclination (see e.g. [Connors fc StarklHOTTl: 



Stark & ConnorsI 119771 IConnors et al. Il980l iLaor et al. 


199C 


; lAgoll 19971: Agol & Krolik 20001: Dovciak et al. 


2004 


. 120081 iLi et all 120091: iSchnittman & KrolikI 120091) . 



Here we briefly review the effect of spin on the polar- 
ization results and refer the reader to previous work for 
further discussion. The solid and dashed curves in Figure 
[in] are computed for a* = and a, = 0.99, respectively. 
Increasing the spin parameter has the well-known effect 
of shifting the spectral peak to higher energies. 
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Fig. 17. — Full disk specific intensity (top panel), degree of po- 
larization (middle panel) and polarization position angle (bottom 
panel), as viewed by an observer at infinity at an inclination of 
73° to the disk normal. All spectra are generated from relativistic 
disk models assuming L/L^^^ = 0.1, a = 0.01, M = IOMq, and 
a* = 0, which are described further in the text. The solid curve 
corresponds to our fiducial model for which B^i^ = B^q/^O = Bq. 
The other curves correspond to different assumption for the mag- 
netic field strength: no magnetic field (triple dot-dashed curve), 
^^ph = Bq/IO (dashed curve), Bpjj = Bo/3 (dotted curve), and 
^^ph = 3i?o (dot-dashed curve). 

The thin curves show the effects of the spacetime on 
polarization most clearly. These models assume the po- 
lari zation at the disk sur face is everywhere identical to 
the lChandrasekha^ (jl960l ) value. Because spacetime near 
the black hole is strongly curved, photons traveling at an 
oblique direction at infinity are in some cases emitted at 
fluid-frame directions rather closer to the polar axis. The 
parallel transport of polarization then leads to an effec- 
tive dilution of the observed polarization. This effect is 
strongest at small radii, where the highest-energy pho- 
tons are predominantly created, so the polarization di- 
minishes somewhat toward higher energy. There is also 
a net rotation of the polarization angle which can be 



seen in the bottom panel. For observer viewing angles 
as shown, which are less than 90° to the disk rotation 
axis (defined in a right-handed sense), the sense of the 
rotation of the polarization position ang le is clockwise 
as on e moves to higher photon energies (IConnors et al.l 
Il980l ). Note, however, that all our calculations omit an- 
other relevant relativistic ef fect: the scattering of return- 
ing radiation. A s shown by|A gol fc Krolik* (200(f) and in 
greater detail bv ISchnittman fc Krolik (2009), this effect 
can lead both to a rise in polarization at high energies 
and a nearly 90° rotation of its direction. 

The thick curves in the bottom two panels of Figure 
116! show the results from disk models for which equation 
([5]) specifies the polarization at the disk surface. The 
low energy polarization is significantly reduced by Fara- 
day depolarization, and is negligible at or below 0.1 keV. 
The degree of polarization increases with photon energy, 
but even near the peak, the polarization is reduced sig- 
nificantly from the non-magnetized, Thomson scattering 
models. 

From Figure [111 it is clear that Faraday rotation has a 
significant effect for the magnetic field strengths adopted 
in our model. However, it may be the case that our as- 
sumed Bph — Beq/'iO = Bq docs not generally hold or 
that our estimate of Bcq from the thin disk model may 
either be higher or lower than in real flows. Therefore, we 
explore the sensitivity to this assumption in Figure [TTI bv 
plotting several spectra with identical spin, mass, accre- 
tion rate, and inclination, but with different assumptions 
for magnetic field strength. 

At the high energy end, all the models asymptote to a 
similar polarization set by the Chandrasekhar value for 
this inclination, but with a slight reduction due to rela- 
tivistic effects. The effects of varying the magnetic field 
are seen at energies near the spectral peak (in vl^) and 
at low energies. Quite generally, the level of polariza- 
tion decreases at all energies as we strengthen the mag- 
netic field. The wavelength dependence is also sensitive 
to field strength. The model with Bph — (triple-dot- 
dashed curve) shows a monotonic decline in polarization 
as photon energy increases. In the models with lower 
but non-zero magnetic field (Sph = Bq/IO, Bo/3), the 
degree of polarization generally increase as photon en- 
ergy increases and reaches a maximum or begins to level 
off at or just below the spectral peak. For even higher 
fields (i?ph = Bq, SBq), the polarization continues rising 
and only levels off out in the Wien tail. 

We can compare our results with those of lGnedin et al.l 
(|2006f ). who also examined the wavelength dependence 
of polarization in accretion disks subject to Faraday ro- 
tation. They computed the polarization from accretion 
disks with purel y vertical magnetic field using the ana- 
lytic models of (jSilant'evI [2002 ). but neglecting the ef- 
fects of relativistic transfer. They reported the long 
wavelength asymptotic dependence of the polarization 
for various assumptions about radial dependence of the 
magnetic field. We find an approximate asymptotic de- 
pendence of p cx with s 0.3. Long wavelength 
photons are predominantly emitted at larger radii where 
relativistic effects are negligible and Pgas > -Prad, so 
Bpb (X Boq oc ( Shaku ra~fc Sunvaevl[Tl73f) . For this 

radial dependence, ICnedin et al.l (|2006l ) find s = 1/3 in 
reasonable agreement with our result. The agreement is 



14 



not surprising since the depolarization parameter (their 
S) is not very sensitive to the f ield geornetry, and the 
asymptotic dependence found bv lSilant'evI (|2002D is very 
similar to that of equation ([5]). 



7. DISCUSSION 



7.1. 



Implications for Spectral Hardening and 
Continuum Spin Estimation 

One of the methods currently being used to try and 
measure the spins of black holes is fitting the ther- 
mal X-ray continuum data to accretion disk models, 
either directly to the models themselves ()Davis et al.l 
I2006t iMiddleton et al.l l2006l or to color-corrected black- 
body disks u sing color correction factors measured frorn 
the models (IShafee et al.l 120061 : iMcChntock et al.l 120061 : 
iLiu et al] l2008f) . All other things being equal, a (pro- 
grade) disk around a more rapidly spinning black hole 
should produce a harder spectrum than around a non- 
rotating back hole. In order for this method to be vi- 
able, the intrinsic hardness of the locally emitted spec- 
trum, which deviates significantly from a blackbody at 
the same effective temperature because of electron scat- 
tering, must be accurately determined from theory. 

All spectral models to date neglect the effects we have 
been exploring in this paper. In particular, magnetic 
vertical support hardens the spectrum and density inho- 
mogeneities soften the spectrum. In the three simulation 
snapshots we examined, these two effects nearly canceled. 
Remarkably, this cancellation occurred even though the 
hardening/softening due to both effects ranged from 

3 — 12% across the different snapshots. This is sugges- 
tive of a more general result, which could arise because 
the magnetic fields providing the support against gravity 
also generate the density inhomogeneities. Nevertheless, 
we caution that it may be still be the case that these 
cancellations were fortuitous, and that in other regimes 
one or the other of these two effects will dominate. 

While we have achieved some success in quantifying 
the uncertainties in the color correction factors of the 
locally emitted spectrum, it should be borne in mind 
that continuum spectral fitting methods to measure black 
hole spin are also subject to u ncertainty in the radia l 
emissivity pr ofile of the di sk. See Beckwith et al. (2008*) , 
IShafee et all (|2008i\ and iNoble et al.. (,2009. ) for recent 
discussions of this issue. 

7.2. Polarization Predictions and Potential Magnetic 
Field Measurements 

The polarization dependence on magnetic field 
strength displayed in Figure ITTI demonstrates that polar- 
ization can, in principle, be used to constrain the mag- 
netic field strength in real accretion flows. 

However, it is important to keep in mind that we have 
neglected several effects. First, we do not include the 
effect of absorption, which will also tend to decrease the 
polarization at the lowest energies, as discussed in Sj5l 
The free-free absorption opacity included in the MC cal- 
culations produces a drop in polarization which occurs 
roughly half a decade lower in energy than the Faraday 
depolarization effect. This is probably an underestimate 
of the effect, since bound-free and bound-bound opacity 
can be significantly greater than free-free opacity for en- 
ergies in the 0.1-2 keV range. The bound- free opacity is 



clearly significant in the 0326c and 0528a equivalent Id 
spectra (Figs. |5]and|6l respectively). However, for the 
hotter annuli which dominate the output near the spec- 
tral peak, the ionization state is sufhciently high that 
bound-free opacity has at most a modest effect (see e.g. 
Fig. HI). 

Faraday rotation and absorption opacity can occasion- 
ally interact in subtle ways to increase the emergent po- 
larization, depending on the vertical gradient of the ther- 
mal source function in the atmosphe re. The phy s ics be - 
hind this is discussed extensively bv lAgol et al.l (|199^ , 
who demonstrated that this can happen near the Balmer 
edge in disk atmosphere models appropriate for quasars. 
It is possible that such effects may also occur near bound- 
free absorption edges in X-ray binary disks. 

Although astronomical X-ray spectropolarimetry has 
hitherto been liniited t o observations of the Crab Nebula 
(jWeisskopf et al.lll976f ). there is now a realistic possibility 
of observing these polarization effects in Galactic black 
holes. The Gravit y and Ext reme Magnetism Small Ex- 
plorer mission ( Ja hoda et al.i 2007: Swa nk et al,.2008" ) is 
now scheduled for launch in 2014 and promises to have 
both the sensitivity and energy resolution in the 2-10 
keV band to detect these effects in a number of objects. 
Similar technologies are also being considered for an X- 
ray polarimetry instrument on board the future Interna- 
tional X-ray Observatory. 

7.3. Generality of these predictions 

If our assumed scaling of _Bph with B^q generally holds, 
Faraday rotation is likely to be important in most lumi- 
nous accretion regimes. The innermost regions of near 
Eddington accretion disks are radiation dominated for 
~ IOMq black holes, and are even more radiation dom- 
inated for supermassive bla ck holes. Using the radia- 
tion dominated relations of IShakura &: SunvaevI ()1973f ) 
we can compute the characteristic rotation angle for op- 
tical depth unity near the spectral peak using equation 

m 



XF,p 



3X1 

-Bph • k. 



(6) 



where Ap ~ ch / (AkBTdf) . The effective temperature 
is computed from the flux 

ic^ThRu 



cff 



Tcs — 



1/4 



(7) 



where Rg = GM/c^ is the gravitational radius, Kcs is 
the electron scattering opacity, a the Stefan-Boltzmann 
constant, r is the radius in gravitational radii, rj is the 
radiative efficiency, m is the accretion rate normalized to 
the Eddington rate MEdd = 47ri?gc/(Kcs'7)j and Rr is a 
function encapsulating both relativistic corrections and 
the ISCO torque boundary condition (notation following 
lKrolikl ll999l). Defining B^r. = (STrPra d)^^^ and using the 
results of IShakura &: SunvaevI ()1973l ). we have 



B, 



B. 
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40 
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\aKcsRgJ 
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K Rr ) 



1/2 



,(8) 



where like summarizes both relativistic and 
boundary condition corrections to the stress profile, and 
i?z is the relativistic correction to the vertical gravity. 
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Combining equations ©-(HI) with the definition of Ap we 
arrive at 



XF,p 



1 



(-) 



1/2 



3/4 
^0 



1/2 2/n 
Q^0.0l'?4<^40 



(9) 



where b is a unit vector in the direction of the mean field, 
q = hc/{kBTcsXp), and Q = B^q/Bpi,. 

Remarkably, the characteristic Faraday rotation angle 
is independent of black hole mass wherever the disk is 
dominated by radiation pressure. Moreover, at least to 
the accuracy of this estimate, it is always ~ 0{1) for 
wavelengths near the thermal peak for the radii where 
most of the light is produced when m ^ 77. Although the 
scaling we have derived should be robust, the actual mag- 
nitude of the effect is subject to significant uncertainty. 
As we have already discussed, the magnitude of a could 
likely be somewhat larger than the fiducial value we have 
chosen for it. Color corrections could alter q. Although 
simulations at a variety of radii have all pointed toward 
Q ~ 40, we cannot say for certain whether this ratio ap- 
plies in the innermost disk. The particular radius dom- 
inating the light output depends on black hole spin and 
the ISCO boundary condition. Because r ~ 10 is appro- 
priate for non-rotating black holes with little torque at 
the ISCO, spin and extra stress in the inner disk would 
diminish the characteristic rotation angle, both by de- 
creasing the characteristic radius of emission and by in- 
creasing i?p{ at small radii (typically Rt — Rn, while 
> !)■ 

We therefore predict that Faraday rotation could sub- 
stantially affect the emergent polarization both for high- 
accretion rate Galactic black holes (in the thermal or 
steep power-law states) and for ma ny AGN (as ha s 
been previously pointed out by, e.g., iBlandford Il990f ) . 
Prospects for using this effect to constrain disk field prop- 
erties are better in the stellar-mass case than in the su- 
permassive case, however, because in AGN it appears 
that the observed polariza tion does not ari se from the 
accretion disk itself (e.g. lAntonuccil 119881 ). but from 
scattering by material further out. Accretion disks in 
cataclysmic variables also emit the bulk of their light 
in the optical and ultraviolet, but these disks are not 
scattering dominated. The theoretically predicted con- 
tinuum polarizatio n in the outburst phase is very small 
([Cheng et al .111 98*81 ) . and has in fact defied observational 
attempts to separat e it out f rom interstellar polarization 
( fNaylor et al. .1996t lMofFat et al. 2001) ■ X-ray polarime- 
try of X-ray binaries might therefore be an optimal way 
of constraining accretion disk magnetic fields. 

8. CONCLUSIONS 

Using snapshots of local shearing box simulations of 
accretion disks with a broad range of radiation to gas 
pressure ratios, we have calculated how a realistic verti- 
cal structure established by MRI turbulence affects the 



emergent photon spectrum and polarization. Most of the 
dissipation within the turbulence occurs at high effec- 
tive optical depth in all the physical regimes that have 
been simulated thus far. As a result, the spectra are 
very insensitive to the details of that dissipation profile, 
and are extremely close to standard model spectra that 
assume that the vertical profile of dissipation per unit 
mass is constant. On the other hand, standard mod- 
els usually neglect the fact that the photospheric regions 
are supported against gravity by magnetic forces. These 
forces reduce the horizontally-averaged densities at the 
thermalization photosphere, resulting in harder Id model 
spectra compared to standard models. 

In addition to lifting the atmosphere, those same 
magnetic forces also produce complex three-dimensional 
density inhomogeneities in the photospheric regions. 
Through a comparison between 3d Monte Carlo calcu- 
lations and Id calculations that were designed to have 
the same emergent fiux, we demonstrated that these in- 
homogeneities act to soften the spectrum. Somewhat 
surprisingly, this softening largely cancels the harden- 
ing due to vertical magnetic support in all three snap- 
shots, even though magnitude of each effect differs be- 
tween the three simulations. Therefore, it may be that 
the color c orrections derived from, e.g. BHSPEC models 
(jPavis fc~H ubcnv 2006), will be approximately correct. 
However, we cannot rule out the possibility that this can- 
cellation was fortuitous in the three snapshots we consid- 
ered here. Given the magnitude of the effects we have 
found, we caution observers that color correction factors 
derived from atmosphere models that neglect magnetic 
support and 3d inhomogeneities may be incorrect by ap- 
proximately ten percent in either direction. 

Because the density inhomogeneities produce a com- 
plex, structured photosphere, the degree of polarization 
of the emerging radiation field is reduced compared to 
a plane-parallel atmosphere, but only slightly. Faraday 
rotation by the strong magnetic fields in the atmosphere 
produces a much stronger effect, producing significant 
reduction in polarization for photon energies near or be- 
low the peak in the local spectrum. This fortunate co- 
incidence that Faraday rotation is strong, but not too 
strong, means that future X-ray polarimetry measure- 
ments of the polarization of the thermal component of 
black hole X-ray binaries could be used as a diagnostic 
of magnetic fields in disks. More extensive theoretical 
calculations will be necessary, however, before this can 
be done quantitatively. 
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APPENDIX 

IMPLEMENTATION OF POLARIZED COMPTON SCATTERING IN THE MONTE CARLO METHOD 



iBerestetskii et ail ()1982[ ) provide an incisive summary and first-principles derivation of the cross sections for both 
polarized and unpolarized Compton scattering. Polarized radiation can be completely described by a Stokes vector 
I = (/, Q, U, V)^ , where / is the total intensity, Q an d U are intensities of linearly polarized radiation, and V is the 
intensity of circular polarization l|Chandrasekhadll960D . We only consider polarization arising from electron scattering. 
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so y = and we will drop it in what follows. Since we only consider scattering of individual photon packets, it will 
be convenient to work in terms of normalized Stokes p arameters q = Q/I and u = U/I. Note that these quantities 
correspond to ^3 and ^1, respectively, in the notation of iBerestetskii et all (|1982D . 

We consider the general problem of modeling electron scattering in the comoving frame of the fluid (as opposed to 
the electron rest frame). Wc define photon and electron momentum four vectors 



i^p =tjmcc(l,k) 

=7771cC(l,/3), 



(Al) 



where to = hv/{meC^) and all other variables have their usual meanings. It is convenient to work with relativistic 
invariants x = 2P^^K^/{mcc)'^ and x' = 2P'^K'^/ {rriec)'^ , where primes denote the scattered quantities. Conservation 
of momentum then requires 



x' — X — 2u!Llj'{1 — cos 6*), 

where cos^ = k • k' is fluid frame scattering angle. 
The differential scattering cross section can be summed over final photon polarizations to yield 



(A2) 



da 



16^ 



- + 4 ) ^ms+m-qn- 

X X ' 



(A3) 



where S = 1/x — 1/x' and gin is the Stokes parameter evaluated in the "internal" basis defined by the scattering 
plane. In general this differs from the "external" basis defined relative to the fluid frame. They are related to each 
other through rotation matrices L{x), which transform I to the internal basis befo re scattering, and £(— yQ, which 
transforms I' back to the external basis after scattering. The matrices are given by (|ChandrasekhaHll960f) 



/I ^ 

Lix) ^ cos2x sin2x 
\0 -sin2x cos2x/ 



(A4) 



U sing these rotation matrices, one flnds qin = q cos 2% — u sin 2% 

iNagirner fc PoutanenI ()1993[ ) provide a thorough discussion of polarization bases and derive expressions for the angles 
X and x' in their appendices. Here, we just quote their results 



cosx = 



smx = 



cos X 



smx 



7(fc^ - cos 9k^)~(3j (1 - cos e)f3^ + k-l3{k'^ - k^) 



7(1 - f3--k){kyk'^ - k^k'y) - /37(1 - COSe)0,ky - Pyk,) 



v/r^yp^^T)7^(l - cos 6*) 
cos6»fc^) ~ (3-1 cos 9) j3^+ k'-^{k^ - k'^] 



(A5) 



y/l^n^y/{6~l)/d{l - cos 9) 

7(1 - f3-h'){k'yk^ - k'^ky) - /37(1 - cos9)0,k'y - pyk'^) 
yT^v/(5~Tp(l - cos 9) 



(A6) 



Here, /? = /3//3 is a unit vector in the electron momentum direction. 

The effects of scattering on the Stokes vector are accounted for by the matrix R = L(x')SL(—x) where F is the trans- 
formation induced by Compton scattering in the internal basis (see e.g. INagirner fc Poutanenl[l993l : IBerestetskii et ahl 
[198^: 



'Sa + SbSa ^ 

Sc Sa 
. Sd, 



(A7) 



where Sa = 4:5(5 + 1) + 2, S'& = [x'/x + x/x') - 2, Sc = Sa - 2, and Sd = 46 + 2. In the ^ hmit, x' 2uj' , 
X 2uj, and 6 (cos 9 — l)/2 so that Sa — >■ cos^ 9 — I, Sb ^ 0, Sc ^ cos^ 9 + 1, Sd ^ 2 cos 9. Therefore, S corresponds 
to the Rayleigh matrix for Thomson scattering in the appropriate limit. Carrying out the matrix multiplication, we 
flnd 



/ Sa + St Sc cos2x -ScSinx ^ 

R= 5cCOs2x' 5*0 cos 2x cos 2x' + S'd sin 2x sin 2x' — 5a sin 2x cos2x' + 'S'd cos 2x sin 2x' 
\~Sc sin 2x' ~Sa cos 2x sin 2x' + 5*^ sin 2x cos 2x' Sa sin 2x sin 2x' + Sd cos 2x cos 2x' / 



(AS) 
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Monte Carlo Implementation 

iPozdniakov et all (|1983D (hereafter PSS83) describe and evaluate Monte Carlo methods for Compton scattering of 
unpolarized radiation. Using the above relations, it is straightforward to generalize the methods presented in PSS83 
to include polarization. For the sake of brevity, we focus only on the additional complexity offered by the polarization 
physics and refer the reader to PSS83 for further details. 

The angle integrated cross section a{x) is independent of polarization, and equivalent to eq. (2.10) of PSS83. 
Therefore, the mean free path and selection of the electron momentum can be evaluated using the exact same method 
as presented in sections 9.4 and 9.5 of PSS83. If the photon momentum is accepted, the scattered photon direction 
k' is then drawn randomly, and cos0 and x' are evaluated as before. However, the scattering probability (eq. 9.8 of 
PSS83) now depends on the Stokes parameters through eq. (jA3|) . The quantity X (defined in §2.1.2 of PSS83) is 
replaced by the quantity in square brackets in eq. (jA3p . Therefore, we must evaluate gin, which requires calculation 
of X via eqs. (lASp . 

If accepted, the photon direction and energy are updated as in PSS83. Additionally, the normalized Stokes param- 
eters after scattering must be calculated. Using I' = R • I and eq. IA8l we find 

, i?2i + i?22'7 + R23U 
q = ^ 

^^11 + ^129 + ^13""' 
, i?3i + i?32(7 + R33U I . , 

u — . (A9j 

^11 + Ri2q + R13U 
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